wineデータはワインの苦味を決める要因に関する実験データである。
考える要因は、ワイン生成時の気温(temp)の2水準(cold or warm)と
果肉と皮の接触の有無(contact)の2水準(no or yes)である。
9人の審査員(judge)が4通りある条件毎に2本ずつ(合計8本)の苦味を評価したデータである。
response : ワインの苦味を表すスコア(0-100点)rating : スコア(response)を5段階評価にグループ分けtemp : ワイン生成時の気温(cold or warm)contact : 果肉と皮の接触の有無(no or yes)bottle : ボトルの識別番号(合計8本)judge : 審査員の識別番号(合計9名)library(ordinal)
## Warning: package 'ordinal' was built under R version 3.0.2
data(wine)
str(wine)
## 'data.frame': 72 obs. of 6 variables:
## $ response: num 36 48 47 67 77 60 83 90 17 22 ...
## $ rating : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 2 3 3 4 4 4 5 5 1 2 ...
## $ temp : Factor w/ 2 levels "cold","warm": 1 1 1 1 2 2 2 2 1 1 ...
## $ contact : Factor w/ 2 levels "no","yes": 1 1 2 2 1 1 2 2 1 1 ...
## $ bottle : Factor w/ 8 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 1 2 ...
## $ judge : Factor w/ 9 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ...
比例オッズ性を仮定した累積ロジットモデル： \[ \log \frac{{\rm Pr}(rating \le j)}{{\rm Pr}(rating > j)} =\theta_j - \beta_1\cdot temp - \beta_2 \cdot contact \]
# 累積ロジットモデル(比例オッズ)
res <- clm(rating ~ temp + contact, data = wine, Hess = TRUE)
summary(res)
## formula: rating ~ temp + contact
## data: wine
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 72 -86.49 184.98 6(0) 4.01e-12 2.7e+01
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## tempwarm 2.503 0.529 4.73 2.2e-06 ***
## contactyes 1.528 0.477 3.21 0.0013 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -1.344 0.517 -2.60
## 2|3 1.251 0.438 2.86
## 3|4 3.467 0.598 5.80
## 4|5 5.006 0.731 6.85
predict関数を用いてnewdataに対する各種予測値行う。
newdata <- unique(wine[c("temp", "contact")])
newdata
## temp contact
## 1 cold no
## 3 cold yes
## 5 warm no
## 7 warm yes
各カテゴリに属する確率：\( {\rm Pr}(rating = j) \), \( j=1,\ldots,5 \)
predict(res, newdata)
## $fit
## 1 2 3 4 5
## 1 0.206790 0.5706 0.1923 0.02362 0.00665
## 2 0.053546 0.3776 0.4431 0.09582 0.02993
## 3 0.020888 0.2014 0.5016 0.20049 0.07563
## 4 0.004608 0.0538 0.3042 0.36360 0.27378
カテゴリの予測値＝属する確率が最大のカテゴリ
predict(res, newdata, type = "class")
## $fit
## [1] 2 3 3 4
## Levels: 1 2 3 4 5
累積確率：\( {\rm Pr}(rating \le j) \), \( j=1,\ldots,5 \)
predict(res, newdata, type = "cum.prob")$cprob1
## 1 2 3 4 5
## 1 0.206790 0.77744 0.9697 0.9933 1
## 2 0.053546 0.43119 0.8743 0.9701 1
## 3 0.020888 0.22230 0.7239 0.9244 1
## 4 0.004608 0.05841 0.3626 0.7262 1
線形予測子：\( \theta_j-\beta_1\cdot temp-\beta_2\cdot contact \), \( j=1,\ldots,5 \)
predict(res, newdata, type = "linear.predictor")$eta1
## 1 2 3 4 5
## 1 -1.344 1.251 3.4669 5.0064 1e+05
## 2 -2.872 -0.277 1.9391 3.4786 1e+05
## 3 -3.847 -1.252 0.9638 2.5033 1e+05
## 4 -5.375 -2.780 -0.5640 0.9755 1e+05
比例オッズ性を仮定した累積ロジット混合モデル：
\[
\log \frac{{\rm Pr}(rating \le j)}{{\rm Pr}(rating > j)}
=\theta_j - \beta_1\cdot temp - \beta_2 \cdot contact + u_{judge},\quad u_{judge} \sim N(0,\sigma^2)
\]
clmm関数の場合はlme4パッケージのlmer関数と同じ方法で柔軟にランダム効果が指定できる。
res <- clmm(rating ~ temp + contact + (1 | judge), data = wine, Hess = TRUE)
summary(res)
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: rating ~ temp + contact + (1 | judge)
## data: wine
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 72 -81.57 177.13 331(996) 1.04e-05 2.8e+01
##
## Random effects:
## Groups Name Variance Std.Dev.
## judge (Intercept) 1.28 1.13
## Number of groups: judge 9
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## tempwarm 3.063 0.595 5.14 2.7e-07 ***
## contactyes 1.835 0.513 3.58 0.00034 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -1.624 0.682 -2.38
## 2|3 1.513 0.604 2.51
## 3|4 4.229 0.809 5.23
## 4|5 6.089 0.972 6.26
clmm2関数の場合はrandom=でランダム切片のみ指定できる。
res <- clmm2(rating ~ temp + contact, random = judge, data = wine, Hess = TRUE)
summary(res)
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## Call:
## clmm2(location = rating ~ temp + contact, random = judge, data = wine,
## Hess = TRUE)
##
## Random effects:
## Var Std.Dev
## judge 1.279 1.131
##
## Location coefficients:
## Estimate Std. Error z value Pr(>|z|)
## tempwarm 3.063 0.595 5.145 2.7e-07
## contactyes 1.835 0.513 3.580 0.00034
##
## No scale coefficients
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 -1.624 0.682 -2.379
## 2|3 1.513 0.604 2.507
## 3|4 4.229 0.809 5.227
## 4|5 6.089 0.972 6.261
##
## log-likelihood: -81.57
## AIC: 177.13
## Condition number of Hessian: 27.62